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ABSTRACT 


We  extend  the  Regional  Seismic  Travel  Time  (RSTT)  tomographic  effort  to  North  America.  In  previous 
work  we  developed  a  real-time  method  to  capture  the  lst-order  effects  of  3 -dimensional  crust  and  upper 
mantle  structure  on  RSTTs.  The  model  parameterization  is  a  global  tessellation  of  nodes  with  a  velocity 
profile  at  each  node.  Interpolation  of  the  velocity  profiles  generates  a  3 -dimensional  crust  and  laterally 
variable  upper  mantle  velocity.  The  upper-mantle  velocity  profile  at  each  node  is  represented  as  a  linear 
velocity  gradient,  which  enables  travel  time  computation  in  approximately  1  millisecond.  Last  computation 
allows  the  model  to  be  used  in  routine  analyses  and  in  operational  monitoring  systems.  Model  velocities  are 
optimized  for  travel-time  prediction  using  a  tomographic  formulation  that  adjusts  the  mantle  velocity  at  the 
Moho,  the  mantle  velocity  gradient,  and  the  average  crustal  velocity.  After  tomography  across  Eurasia, 
rigorous  tests  find  that  Pn  travel  time  residuals  are  reduced  from  a  standard  deviation  of  approximately 
1.75  seconds  (akl35  model)  to  approximately  1.25  seconds,  further,  location  error  is  consistently  reduced 
by  approximately  45%  for  events  located  using  the  Pn  phase.  Lor  North  American  tomography  we  begin  by 
reconciling  North  American  seismic  bulletins.  To  these  bulletins  we  add  high-quality  data  sets  that  improve 
data  coverage  and  data  quality,  including  arrival  times  reported  by  the  Array  Network  facility  for  US  Array 
stations.  In  the  western  United  States,  US  Array  provides  unprecedented  data  coverage  with  station  spacing 
of  approximately  70  km.  The  National  Earthquake  Information  Center  (NEIC)  researchers  contribute 
unique  ground-truth  data  sets  that  have  been  culled  over  decades.  NEIC  data  is  particularly  helpful  in 
eastern  North  America,  where  the  seismicity  rate  is  low.  Arrival-time  measurements  from  all  sources  are 
re-associated,  and  event  locations  that  are  best  constrained  by  seismic  data  are  relocated  using  all  available 
arrivals.  Events  with  known  locations,  e.g.,  explosions,  are  not  relocated  and  these  events  help  to  establish 
absolute  travel  time  accuracy.  Seismic  locations  are  evaluated  against  network  coverage  criteria  to  estimate 
hypocenter  accuracy.  The  full  error  budget  for  each  travel  time  -  hypocenter  uncertainty  and  arrival-time 
measurement  uncertainty  -  is  evaluated  to  provide  a  datum-specific  uncertainty  that  establishes  weighting 
in  the  tomographic  inversion.  Model  validation  includes  prediction  of  travel  times,  as  well  as  relocation  of 
ground-truth  events  to  measure  location  accuracy.  Neither  the  travel  time  validation  data  nor  the  location 
validation  data  are  used  in  the  tomographic  inversion. 
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OBJECTIVES 

Monitoring  seismic  events  at  ever-lower  magnitude  thresholds  requires  the  utilization  of  seismic  stations 
that  are  close  to  the  event.  The  proliferation  of  seismic  stations  across  the  globe  results  in  lower  seismic 
detection  thresholds  because  more  stations  are  likely  to  be  at  regional  distance  to  an  event.  However, 
utilization  of  arrival  time  data  at  regional  distances  often  results  in  degraded  location  accuracy,  because 
large  variations  in  regional  crust  and  upper-mantle  structure  can  result  in  large  travel  time  prediction  errors. 

In  previous  work  we  successfully  developed  an  operations-oriented  method  to  account  for  the  1  st-order 
effects  of  3 -dimensional  crust  and  upper  mantle  structure  on  RSTTs  (Myers  et  al.,  2010).  The  RSTT  model 
parameterization  is  a  global  tessellation  of  nodes,  so  any  model  in  RSTT  format  is  inherently  global  in 
extent  and  a  regional  travel  time  can  be  computed  anywhere  for  the  dominant  regional  phases  ( Pn ,  Pg,  Sn , 
and  Lg).  However,  in  order  to  improve  travel  time  prediction  accuracy  over  conventional  methods,  the 
model  must  be  “tuned”  along  each  ray  path.  RSTT  model  tuning  is  accomplished  by  first  constructing  a 
best  approximation  to  the  3-dimensional  structure  of  the  Earth’s  crust  and  upper-most  mantle  based  on 
existing  models.  Model  velocity  -  but  not  layer  interfaces  -  are  then  adjusted  using  tomographic  imaging 
methods  to  further  improve  prediction  accuracy.  RSTT  model  refinement  was  conducted  first  over  the  -1/4 
of  the  globe  containing  Eurasia  and  North  Africa  (Myers  et  ah,  2010;  Figure  1).  RSTT  tomography  across 
Eurasia  and  North  Africa  results  in  a  factor  -2  improvement  in  location  accuracy. 

Here  we  extend  RSTT  tomography  to  North  America.  By  extending  tomography  to  North  America  we 
address  the  need  to  improve  global  tomographic  coverage,  which  is  desirable  for  CTBT  monitoring,  and  we 
motivate  RSTT  use  at  the  NEIC,  whose  primary  interest  is  monitoring  earthquakes  in  the  continental  U.S. 


Figure  1.  Phased  approach  for  applying  RSTT  tomography,  with  the  end  goal  of  producing  a  global 
model  for  the  universal  prediction  of  regional  phases.  Phase  1  is  complete,  and  Phase  2 
would  be  initiated  under  this  proposal.  (Reproduced  from  Myers  et  al.,  2010) 
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RESEARCH  ACCOMPLISHED 

Work  Plan 

RSTT  model  refinement  has  5  essential  components.  First,  based  on  previously  determined  models,  a  prior 
model  is  built  that  includes  lateral  variability  of  the  geologic  layers  that  comprise  the  crust.  The  prior  model 
also  includes  lateral  variability  of  P- wave  and  S-wave  velocity  in  the  crustal  layers  and  in  the  upper  mantle. 
A  good  prior  model  is  essential  because  the  RSTT  tomography  inversion  is  invariably  under  determined, 
necessitating  a  model  search  in  the  neighborhood  of  the  starting  model.  Second,  a  high-quality  data  set  of 
event  locations  and  arrival-time  observations  is  compiled  for  use  in  tomography  and  for  model  validation. 
Third,  tomographic  inversion  is  used  to  adjust  model  velocities  so  that  travel  time  predictions  are  in 
agreement  with  the  high-quality  data  set.  Layer  boundaries  are  not  adjusted  in  the  tomographic  inversion 
because  RSTT  data  sets  cannot  adequately  resolve  the  trade-off  between  boundary  depth  and  velocity. 
Fourth,  assessment  of  travel  time  error  is  accomplished  by  comparing  model  predictions  to  the  validation 
data  set  (not  used  in  tomography).  The  fifth  step  is  model  validation,  in  which  events  with  known  locations 
are  relocated  to  measure  location  accuracy  and  compare  observed  location  error  with  estimates  derived 
from  model  uncertainty. 

To  date,  we  have  completed  the  data  compilation  effort,  begun  evaluation  of  the  CRUST2.0  starting  model, 
and  conducted  preliminary  tomographic  inversions.  We  report  on  the  data  set  and  on  our  preliminary 
findings. 

Model  Parameterization 

Myers  et  al.  (2010)  provide  a  detailed  description  of  RSTT  model  parameterization,  travel  time  calculation, 
and  tomographic  formulation,  and  we  provide  a  brief  review  here.  Crust  and  upper  mantle  velocity 
structure  are  represented  using  radial  velocity  profiles  at  geographically  distributed  nodes  (Figure  2).  The 
nodes  form  a  triangular  tessellation  that  seamlessly  covers  the  globe.  Node  spacing  is  approximately  1°  for 
the  model  presented  here.  Node  spacing  may  be  adjusted  as  needed,  and  we  plan  to  investigate  the  use  of 
0.25°  node  spacing  in  North  America  that  may  be  warranted  by  outstanding  data  coverage  in  the  western 
U.S.  Velocity  interfaces  are  defined  by  the  radial  distance  from  the  center  of  the  Earth,  which  allows  us  to 
explicitly  build  the  GRS80  ellipsoid  (Moritz,  1980)  into  the  model  and  obviate  travel-time  corrections  for 
ellipticity. 

We  adopt  the  velocity  versus  depth  profile  in  the  crust  from  Pasyanos  et  al.  (2004),  which  includes  model 
layers  for  water,  3  types  of  sediments,  upper  crystalline  crust,  middle  crust,  and  lower  crust  (Figure  2).  The 
crustal  layers  overlay  a  mantle  velocity  profile  that  is  simplified  to  two  parameters:  velocity  at  the  Moho 
and  a  linear  velocity  gradient  with  depth.  By  interpolating  model  parameters  from  surrounding 
nodes  -  layer  thickness,  velocity,  and  mantle  gradient  -  we  generate  a  continuous  model  of  the  3-D  crust 
and  laterally  varying  upper  mantle. 


Travel  Time  Calculation 

Parameterization  of  upper  mantle  velocity  with  a  linear  gradient  facilitates  an  approximation  for  Pn  travel 
time  that  enables  real-time  computation  (-1  millisecond).  Computation  of  Pn  travel  time  at  near-regional 
distance  (<700  km)  commonly  assumes  that  the  Pn  phase  propagates  as  a  head  wave,  with  a  ray-path  that 
follows  the  contour  of  the  Moho  (e.g.  Hearn,  1984).  The  head  wave  assumption  results  in  poor  travel  time 
prediction  at  far-regional  distance  (>-700  km)  because  the  Pn  ray  can  dive  appreciably  into  the  mantle  due 
to  a  positive  velocity  gradient  with  depth  and  Earth  sphericity  (e.g,.  Zhao  and  Xie,  1993;  Ritzwoller  et  al., 
2003;  Heam  et  al.,  2004).  To  more  accurately  predict  Pn  at  far-regional  distances,  Zhao  (1993)  and  Zhao 
and  Xie  (1993)  employ  a  constant  linear  velocity  gradient  in  the  upper  mantle  for  the  whole  study  area. 
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a)  Geographic  node  sampling 


b)  Model  parameterization 
at  each  node 
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Figure  2.  Global  model  parameterization,  a)  An  example  tessellation  with  approximately  5°  grid 
spacing.  The  inset  shows  the  1°  used  in  this  study.  Color  indicates  Moho  depth  of  the 
starting  model,  b)  An  example  velocity/depth  profile  as  defined  at  each  node.  The  mantle 
portion  of  the  profile  is  specified  by  the  velocity  at  the  crust/mantle  interface  and  a  linear 
gradient.  (Reproduced  from  Myers  et  al.,  2010) 

The  Zhao  (1993)  and  Zhao  and  Xie  (1993)  travel-time  calculation  is  similar  to  the  widely  used  approach  of 
Hearn  (1984),  with  an  additional  term  (y)  introduced  to  account  for  diving  rays.  The  travel-time  calculation 
is 


TT  =  idiSi+a  +  j3  +  r  (1) 

i=  1 

where  d  and  s  are  the  distance  and  slowness  (taken  as  1 /velocity  below  the  Moho)  in  each  of  the  i  segments 
comprising  the  great-circle  path  between  Moho  pierce  points  near  the  event  and  station,  a  and  p  are  the 
crustal  travel  times  at  the  source  and  receiver,  and  y  (described  below). 


From  Zhao  (1993)  and  Zhao  and  Xie  (1993), 


r= 


C2Xl 

24  Vo 


(2) 


where  Xm  is  the  horizontal  distance  traveled  in  the  mantle,  and  V0  is  a  regional  average  of  mantle  velocity  at 
the  Moho.  c  =  g*  s  +l/r,  where  1/r  is  an  Earth  flattening  correction  and  r  is  the  radius  at  which  a  ray 
enters  and  exits  the  linear  velocity  gradient,  g  (Helmberger,  1973;  Zhao  and  Xie,  1993).  This 
approximation  is  valid  when  ch«  1,  where  h  is  the  bottoming  depth  of  the  ray  in  a  linear  velocity  gradient. 


We  use  a  spatially  varying  mantle  velocity  gradient,  c  (Phillips  et  al.,  2007),  and  we  calculate  y  by 
averaging  c  along  the  ray  track.  V0  remains  an  average  Pn  velocity  over  the  whole  model,  which  allows  us 
to  take  advantage  of  linear  tomographic  inversion  methods  (see  below).  Tests  find  that  using  a  global 
average  for  V0  introduces  negligible  travel-time  error  when  Pn  velocities  range  from  7.5  km/s  to  8.3  km/s. 


Data 

Ongoing  GNEM  R&D  integration  efforts  include  the  collection  of  seismic  bulletins  that  comprise  the  vast 
majority  of  events  and  arrival  time  measurements  that  are  available  in  North  America.  To  these  ongoing 
efforts  we  add  additional  high-quality  data  sets  that  provide  coverage  in  specific  areas.  The  USarray  project 
(http ://www. iris . edu/U S Array/researchers/data.html)  is  providing  unprecedented  coverage  with  station 
spacing  of  approximately  70  km.  The  transportable  USarray  is  a  staged  deployment  that  began  on  the  west 
coast  of  the  U.S.  and  is  steadily  moving  eastward.  USarray  currently  straddles  the  Rocky  Mountain  front 
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and  the  Great  Plains.  The  transportable  array  is  scheduled  to  reach  the  Eastern  US  in  2013.  Because  of  the 
dense  station  spacing,  arrival  time  measurements  of  regional  distance  events  at  US  Array  stations  can 
constrain  short  wavelength  lateral  variations  in  the  crust  and  upper  mantle.  We  have  loaded  the  arrival  time 
measurements  from  the  Array  Network  Facility  (http://anf.ucsd.edu)  into  the  GNEM  R&D  database  and 
reconcile  the  events/arrivals  with  existing  holdings.  Second,  NEIC  researchers  have  collected  unique 
ground  truth  data  sets  in  North  America,  particularly  in  the  Eastern  U.S.  The  NEIC  data  have  been 
integrated/reconciled  with  GNEM  R&D  holdings  in  a  similar  fashion  as  the  ANF  data. 

Data  processing  follows  the  procedure  outlined  in  Myers  et  al.  (2010).  Epicenter  accuracy  for  each  event  in 
the  reconciled  bulletin  is  assessed  using  the  network  coverage  criteria  of  Bondar  et  al.  (2004).  We  further 
add  non-seismic  constraints  based  on  known  explosion  locations,  ground  displacement  from  interferometric 
synthetic  aperture  radar  (InSAR),  and  mine  locations.  To  diminish  the  possibility  of  introducing  travel 
times  for  phases  that  interact  with  velocity  discontinuities  at  -410  km  and  -660  km,  the  maximum  event- 
station  distance  range  is  set  to  15°.  The  minimum  event-station  distance  range  is  determined  by  the  post- 
critical  refraction  for  a  wave  interacting  with  the  Moho.  In  practice,  the  minimum  distance  varies  from  tens 
of  km  in  the  ocean  (thin  crust)  to  approximately  200  km  in  regions  of  thick  crust. 

Because  the  goal  of  this  work  is  to  produce  a  model  for  Pn  travel-time  prediction  for  real-time  monitoring, 
it  is  important  that  Pn  prediction  error  is  unbiased  relative  to  teleseismic  P- wave  prediction  error.  Previous 
efforts  have  achieved  unbiased  Pn  error  by  using  an  ad  hoc  travel  time  correction  (Yang  et  al.,  2004).  To 
achieve  unbiased  Pn  error,  we  recomputed  each  event  origin  time  using  at  least  10  P- wave  arrivals.  The 
hypocenter  is  then  fixed  during  the  tomographic  procedure,  which  forces  Pn  prediction  error  to  be  unbiased 
relative  to  teleseismic  P-wave  error. 

All  picks  are  evaluated  against  an  uncertainty  budget  that  accounts  for  event  mislocation,  a  global  average 
of  ak!35  prediction  uncertainty,  and  arrival-time  measurement  uncertainty  (See  Myers  et  al.,  2010). 
Following  removal  of  local  outliers,  we  produce  summary  rays  for  each  station  by  averaging  residuals  for 
events  that  are  within  0.25°  of  one  another  (Figure  3). 

For  North  America  we  begin  with  80728  Pn  paths.  Outlier  removal  reduces  the  number  of  paths  to  705 14, 
of  which  7575  are  set  aside  for  validation.  Use  of  summary  rays  reduces  the  number  of  paths  from  62939  to 
35736.  The  standard  deviation  of  the  initial  Pn  data  set  (80728  paths)  is  2.71  seconds  with  a  mean  of  0.81 
seconds.  After  outlier  remove  the  standard  deviation  is  reduced  to  2.35  seconds  with  a  mean  of  0.78 
seconds. 


Figure  3.  RSTT  Pn  summary  ray  paths  for  North  America.  Stations  are  red  triangles  and  events  are 
yellow  “pluses”. 
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Tomography 

The  Pn  travel  time  (Equation  1)  lends  itself  to  a  linear  tomographic  formulation.  Because  our  primary 
objective  is  to  improve  travel-time  prediction,  we  avoid  the  use  of  parameters  that  would  not  be  part  of  a 
subsequent  travel-time  calculation  (e.g.  event  and  station  time  terms).  In  matrix  form,  the  tomographic 
system  of  equations  is: 
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where 

t  =  travel  time 

s  =  mantle  slowness  below  the  Moho  (a.k.a.  Pn  slowness) 

x  =  Pn  distance  (or  weight)  for  each  model  node 

c  =  normalized  velocity  gradient,  v=vo(l  +cz ) 

Xm  =  length  of  Pn  ray  path  in  the  mantle 
Vo  =  average  Pn  velocity 

v  =  velocity  of  a  crustal  layer 

k  =  index  on  K paths  (travel-time  observations) 

p  =  index  on  Q  crustal  layers 

l  =  length  of  the  ray  path  in  a  specified  crustal  layer  (determined  by  layer  thickness  and  ray 

parameter  in  Equations  (2)  and  (3)). 

a  =  node-specific  adjustment  to  the  slowness  of  each  crustal  layer  (crustal  modifier). 


The  tomographic  equation  solves  for  the  model  slowness  below  the  Moho,  s  (a.k.a.  Pn  slowness),  the 
square  of  mantle  velocity  gradient,  c2,  and  a  scalar  adjustment  to  crustal  slowness,  a.  The  formulation  in 
Equation  3  is  similar  to  the  approach  presented  in  Phillips  et  al.  (2007),  with  the  difference  that  we  use  a 
scalar  adjustment  to  the  slowness  of  the  crustal  stack,  as  opposed  to  a  time  term,  to  account  for  travel-time 
errors  in  crustal  legs  of  the  Pn  ray.  The  crustal  legs  can  impart  as  much  or  more  error  on  the  travel  time 
prediction  as  the  travel  time  in  the  mantle.  Because  our  goal  is  to  accurately  predict  subsequent  travel  times 
using  the  tomographic  model,  it  is  important  to  fold  all  adjustments  affecting  travel  time  into  the  model, 
rather  than  absorb  the  error  in  a  time  term  that  is  discarded  and  will  not  be  used  in  subsequent  travel  time 
prediction. 

Preliminary  Results 

Initial  tomgraphic  inversions  for  North  America  are  consistent  with  the  broad  trends  in  previous  studies 
(Figures  4,  5,  and  6).  The  mantle  velocity  at  the  Moho  approaches  8.3  km/s  across  much  of  the  Canadian 
Shield  and  drops  to  between  8.0  km/s  and  8.15  km/s  in  the  Eastern  U.S.  (Figure  4).  West  of  the  Rocky 
Mountain  Front,  velocity  drops  abruptly  below  8  km/s,  with  low-velocity  anomalies  of  approximately 
7.7  km/s  in  some  places. 
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Figure  4.  Mantle  velocity  at  the  Moho.  RSTT  tomography  is  consistent  with  the  known  pattern  of 
upper-mantle  velocity:  very  fast  across  the  Canadian  Shield;  fast  across  the  eastern  U.S., 
and  slow  west  of  the  Rocky  Mountain  Front. 

The  pattern  of  velocity  gradient  in  the  mantle  is  dominated  by  oceanic  subduction  (Figure  5).  Subduction  of 
high  velocity  oceanic  lithosphere  creates  complex  three-dimensional  structure  that  cannot  be  fully  captured 
by  the  RSTT  model  parameterization.  However,  the  first-order  effect  on  Pn  travel  times  is  early  arrival  at 
far-regional  distances  when  rays  dive  into  the  mantle  and  preferentially  channel  through  the  high-velocity 
ocean  lithosphere.  RSTT  captures  this  travel-time  effect  with  a  strong  velocity  gradient,  as  seen  in  southern 
Alaska,  the  Pacific  Northwest,  southern  Mexico,  and  the  Caribbean.  Other  areas  of  high  gradient  include 
the  Wyoming  Archean  terrain  and  the  central  U.S.  We  expect  the  Canadian  Shield  have  a  strong  velocity 
gradient,  but  the  tomographic  data  set,  which  only  has  a  few  long  paths  across  the  Canadian  Shield,  is  not 
suited  to  resolve  mantle  gradient  in  the  Canadian  Shield.  In  subsequent  inversions  we  may  force  a  high 
gradient  in  the  Canadian  Shield. 
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Figure  5.  Mantle  velocity  gradient.  Subduction  zones  are  characterized  with  a  high  velocity  gradient. 

Figure  6  shows  the  modification  to  crustal  slowness  {a  in  equation  3).  The  Pn  data  set  does  not  provide 
resolution  of  each  of  the  7  layers  that  comprise  the  RSTT  crust,  so  we  solve  for  a  single  modification  to 
crustal  slowness  at  each  node  that  is  uniformly  applied  each  crustal  layer.  Increases  in  crustal  velocity 
(negative  a  values)  are  seen  in  California  and  Idaho.  Decreases  in  crustal  velocity  are  seen  in  the  Colorado 
Plateau,  Wyoming,  and  the  Wrangle  Mountains  of  Alaska.  The  crustal  velocity  does  not  change 
significantly  throughout  eastern  North  America. 

1 90 '  200*  21 0‘  220'  230' 240  260  260' 270'  280 ‘  290' 300*  31  O' 


0.85  0.97  0.89  1.01  1.03  1.05 

Modifier 


Figure  6.  Modification  of  crustal  slowness.  This  term  adjusts  the  slowness  (1/velocity)  of  the  whole 
crustal  stack. 
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Validation 

We  leave  out  10%  of  the  tomographic  data  for  use  in  non-circular  validation  tests.  The  validation  data 
provides  sampling  across  North  America,  so  residual  summary  statistics  are  a  good  measure  of  expected 
model  performance  in  monitoring  systems.  Figure  7  shows  travel  time  error  as  a  function  of  distance  for  the 
akl35  model,  and  the  RSTT  model.  In  this  case  error  is  represented  as  the  average  residual  in  successive 
distance  bins.  Similar  to  the  Eurasia  results,  akl35  uncertainty  is  approximately  2  seconds  at  near-regional 
distance  and  increases  rapidly  beyond  5°  to  a  maximum  of  almost  3  seconds  at  far-regional  distance.  The 
RSTT  model  is  significantly  improved  over  akl35  by  reducing  the  overall  travel  time  prediction  error  down 
to  approximately  1.5  seconds.  The  trend  of  increasing  residuals  as  a  function  of  distance  is  still  evident  in 
the  preliminary  model. 


Figure  7.  Pn  travel  time  prediction  error  as  a  function  of  distance  for  North  American  validation 
data  set:  black  line  is  ak!35,  red  line  is  RSTT. 


CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  completed  the  data  compilation  needed  to  extend  RSTT  tomography  to  North  America.  That 
effort  includes  integration  into  our  database  of  data  provided  by  the  NEIC  and  picks  provided  by  the  Array 
Network  Facility  for  US  Array  stations.  The  integrated  data  have  been  relocated  and  epicenter  accuracy 
criteria  applied.  Pick  outliers  have  been  identified  based  on  overall  statistics  and  comparison  with 
neighboring  data.  Summary  rays  for  the  Pn  data  (Figure  3)  provide  excellent  ray  coverage  across  western 
North  America,  good  coverage  across  the  eastern  U.S.,  poor  coverage  across  central  and  eastern  Canada. 

Preliminary  tomographic  results  use  a  starting  model  that  is  based  on  CRUST2.0.  We  find  that  modification 
to  CRUST2.0  are  needed.  Nonetheless,  preliminary  tomographic  results  are  consistent  with  known 
velocities  for  major  tectonic  provinces.  Tests  of  travel-time  prediction  based  on  a  validation  data  set 
(Figure7)  show  approximately  25%  reduction  in  the  standard  deviation  of  residuals. 
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